数据包加载

knitr::opts_chunk$set(echo = TRUE)
# 检测和安装依赖包
package_list <- c("randomForest","ggplot2","pheatmap")
# 判断R包加载是否成功来决定是否安装后再加载
for(p in package_list){
  if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
    install.packages(p)
    suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
  }
}
library(pheatmap)

文件输入+自定义分组

# 1. 读取OTU表
otutab = read.table("C:/16S/result/randomforestclassification/otutab.txt", header=T, row.names= 1, sep="\t", comment.char = "", stringsAsFactors = F)
# 2. 读取物种注释
tax = read.table("C:/16S/result/randomforestclassification/taxonomy.txt", header=T, row.names= 1, sep="\t",comment.char = "", stringsAsFactors = F) 
# 数据按实验设计手动筛选
metadata_all = read.table("C:/16S/result/randomforestclassification/metadata.txt", header=T, row.names=1, sep="\t") 

SA = ""
SB = ""
SC = ""
Site_A = "S"
Site_B = "E"
Site_C = "G"
State_A = "H"
State_B = "D"
State_C = "I"
Sday_A = "pre"
Sday_B = "med"
Sday_C = "pre"

filetial1 = "rf.txt"
filetial2 = "rf.pdf"
sub_metadatanames = "sub_metadata"
rfcv_names = "rfcv"
importance_names = "importance"
top_f = "top_feautres"
prediction_name = "prediction_binary"
test_names = "test"
class_names = "new"
result_name = "result"
output_site = "C:/16S/result/randomforestclassification/"

result_file_names <- paste(result_name,Site_A,Site_B,"new",sep ="_")
result_file_names <- paste(output_site,result_file_names,sep ="")

sub_metadata_filename <- paste(sub_metadatanames,Site_A,Site_B,filetial1,sep ="_")
sub_metadata_filename <- paste(output_site,sub_metadata_filename,sep ="")

rfcv_filename_txt <- paste(class_names,rfcv_names,Site_A,Site_B,filetial1,sep ="_")
rfcv_filename_txt <- paste(output_site,rfcv_filename_txt,sep ="")

rfcv_filename_pdf <- paste(class_names,rfcv_names,Site_A,Site_B,filetial2,sep ="_")
rfcv_filename_pdf <- paste(output_site,rfcv_filename_pdf,sep ="")

importance_filename_txt <- paste(class_names,importance_names,Site_A,Site_B,filetial1,sep ="_")
importance_filename_txt <- paste(output_site,importance_filename_txt,sep ="")

top_feautre_filename <- paste(class_names,top_f,Site_A,Site_B,filetial2,sep ="_")
top_feautre_filename <- paste(output_site,top_feautre_filename,sep ="")

prediction_filename <- paste(class_names,prediction_name,Site_A,Site_B,filetial1,sep ="_")
prediction_filename <- paste(output_site,prediction_filename,sep ="")

testfilename1 <- paste(class_names,test_names,Site_A,filetial2,sep ="_")
testfilename1 <- paste(output_site,testfilename1,sep ="")

testfilename2 <- paste(class_names,test_names,Site_B,filetial2,sep ="_")
testfilename2 <- paste(output_site,testfilename2,sep ="")

print(result_file_names)
## [1] "C:/16S/result/randomforestclassification/result_S_E_new"
print(sub_metadata_filename)
## [1] "C:/16S/result/randomforestclassification/sub_metadata_S_E_rf.txt"
print(rfcv_filename_txt)
## [1] "C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.txt"
print(rfcv_filename_pdf)
## [1] "C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.pdf"
print(importance_filename_txt)
## [1] "C:/16S/result/randomforestclassification/new_importance_S_E_rf.txt"
print(top_feautre_filename)
## [1] "C:/16S/result/randomforestclassification/new_top_feautres_S_E_rf.pdf"
print(prediction_filename)
## [1] "C:/16S/result/randomforestclassification/new_prediction_binary_S_E_rf.txt"
print(testfilename1)
## [1] "C:/16S/result/randomforestclassification/new_test_S_rf.pdf"
print(testfilename2)
## [1] "C:/16S/result/randomforestclassification/new_test_E_rf.pdf"
# 筛选自己想要的分组以及内容
# sub_metadata1 = subset(metadata_all, State %in% c(State_A))
# sub_metadata1
sub_metadata2 = subset(metadata_all, Site_1 %in% c(Site_A,Site_B))
sub_metadata2
# 实验设计与输入文件交叉筛选。
idx = rownames(sub_metadata2) %in% colnames(otutab)
sub_metadata3 = sub_metadata2[idx,]
sub_metadata3
sub_otutab_tax = otutab[,rownames(sub_metadata3)]
sub_otutab_tax
write.table(sub_metadata3,file = sub_metadata_filename,quote = F,sep = '\t', row.names = T, col.names = T)

# OTU丰度筛选阈值,默认0.1%,0为来筛选
thre = 0.02
# 输出文件名前缀
prefix = "tax_S_E_new"

# 生成各分类级汇总特征表
suppressWarnings(suppressMessages(library(amplicon)))
format2stamp(sub_otutab_tax, tax, thre, prefix)
## Warning: 程辑包'dplyr'是用R版本4.3.2 来建造的
now_dir <- getwd()  # 当前工作目录
aim_dir <- "C:/16S/result/randomforestclassification"

# 创建目标文件夹(如果不存在)
if (!dir.exists(aim_dir)) {
  dir.create(aim_dir, recursive = TRUE)
}

file_patterns <- c("tax_S_E_new_1Kingdom", "tax_S_E_new_2Phylum", "tax_S_E_new_3Class", "tax_S_E_new_4Order", 
                   "tax_S_E_new_5Family", "tax_S_E_new_6Genus", "tax_S_E_new_7Species", "tax_S_E_new_8OTU0.02")

# 遍历文件模式并移动文件
sapply(file_patterns, function(pattern) {
  files_to_move <- list.files(now_dir, pattern = pattern, full.names = TRUE)
  sapply(files_to_move, function(file) {
    file.copy(file, aim_dir, overwrite = TRUE)
    file.remove(file)
  })
})
## tax_S_E_new_1Kingdom.C:/16S/All_figures_pipeline/tax_S_E_new_1Kingdom.txt 
##                                                                      TRUE 
##   tax_S_E_new_2Phylum.C:/16S/All_figures_pipeline/tax_S_E_new_2Phylum.txt 
##                                                                      TRUE 
##     tax_S_E_new_3Class.C:/16S/All_figures_pipeline/tax_S_E_new_3Class.txt 
##                                                                      TRUE 
##     tax_S_E_new_4Order.C:/16S/All_figures_pipeline/tax_S_E_new_4Order.txt 
##                                                                      TRUE 
##   tax_S_E_new_5Family.C:/16S/All_figures_pipeline/tax_S_E_new_5Family.txt 
##                                                                      TRUE 
##     tax_S_E_new_6Genus.C:/16S/All_figures_pipeline/tax_S_E_new_6Genus.txt 
##                                                                      TRUE 
## tax_S_E_new_7Species.C:/16S/All_figures_pipeline/tax_S_E_new_7Species.txt 
##                                                                      TRUE 
## tax_S_E_new_8OTU0.02.C:/16S/All_figures_pipeline/tax_S_E_new_8OTU0.02.txt 
##                                                                      TRUE
# 打印成功信息
print("Files generated and moved to the target directory successfully.")
## [1] "Files generated and moved to the target directory successfully."
# 在指定目录生成tax_1-8共7个级别+OTU过滤文件

分类级选择

先使用format2stamp.Rmd基于OTU表(otutab.txt)、物种注释(taxonomy.txt)和元数据(metadata.txt)筛选样本、高丰度特征,并分类汇总各分类级(tax_1-8)。 然后对各分类级进行准确性评估

# 读取实验设计、和物种分类文件
# 实验用的数据重命名
metadata_f = sub_metadata3
metadata_f
# 设置训练集
train = subset(metadata_f, Site_1 %in% c("S","E"))
train
idx3 = sample(1:386, size = 200, replace = F)
# 注意sample后面的比例要根据上面文件自己调。size就是比例
# 选择的样本标记为TRUE,未选择的为FALSE
idx3 = 1:386 %in% idx3
#train2 = metadata_f[idx3,]
train2 = train[idx3,] # train2是根据实验设计在train里面又筛选的作为训练集的数据。测试集是剩余部分。
train2
train2$Site_1 = as.factor(train2$Site_1)
summary(train2)
##    Filename          Group_day          Group_day1         Group_dayE       
##  Length:200         Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   Group_dayS           Group            StateSite         StateSite_1       
##  Length:200         Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##     State               Site           Site_1       Day        
##  Length:200         Length:200         E: 73   Min.   : 0.000  
##  Class :character   Class :character   S:127   1st Qu.: 0.000  
##  Mode  :character   Mode  :character           Median : 3.000  
##                                                Mean   : 3.435  
##                                                3rd Qu.: 7.000  
##                                                Max.   :10.000  
##     Day_num           Day1               DayE            Individual       
##  Min.   : 0.000   Length:200         Length:200         Length:200        
##  1st Qu.: 0.000   Class :character   Class :character   Class :character  
##  Median : 3.000   Mode  :character   Mode  :character   Mode  :character  
##  Mean   : 3.435                                                           
##  3rd Qu.: 7.000                                                           
##  Max.   :10.000
# 筛选一个合适的分类等级
library(randomForest)
# 自己运行
for(i in c("2Phylum","3Class","4Order","5Family","6Genus","8OTU0.02")){
  # i="5Family"
  set.seed(0)
  table_f = read.table(paste0("C:/16S/result/randomforestclassification/","tax_S_E_new_",i,".txt"),header = T, row.names = 1)
  table_f = table_f[,rownames(train2)]
  train2$Site_1 = as.factor(train2$Site_1)
  rf = randomForest(t(table_f), train2$Site_1, importance=T, proximity=T, ntree = 5000)
  print(i)
  print(rf)
}
## [1] "2Phylum"
## 
## Call:
##  randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000,      importance = T, proximity = T) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 30%
## Confusion matrix:
##    E   S class.error
## E 29  44   0.6027397
## S 16 111   0.1259843
## [1] "3Class"
## 
## Call:
##  randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000,      importance = T, proximity = T) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 29.5%
## Confusion matrix:
##    E   S class.error
## E 30  43   0.5890411
## S 16 111   0.1259843
## [1] "4Order"
## 
## Call:
##  randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000,      importance = T, proximity = T) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 5
## 
##         OOB estimate of  error rate: 22.5%
## Confusion matrix:
##    E   S class.error
## E 34  39  0.53424658
## S  6 121  0.04724409
## [1] "5Family"
## 
## Call:
##  randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000,      importance = T, proximity = T) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 7
## 
##         OOB estimate of  error rate: 23.5%
## Confusion matrix:
##    E   S class.error
## E 33  40  0.54794521
## S  7 120  0.05511811
## [1] "6Genus"
## 
## Call:
##  randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000,      importance = T, proximity = T) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 18%
## Confusion matrix:
##    E   S class.error
## E 39  34  0.46575342
## S  2 125  0.01574803
## [1] "8OTU0.02"
## 
## Call:
##  randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000,      importance = T, proximity = T) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 23
## 
##         OOB estimate of  error rate: 13.5%
## Confusion matrix:
##    E   S class.error
## E 47  26 0.356164384
## S  1 126 0.007874016
# OOB estimate of  error rate: 越小越好。

最佳水平数据读取和统计

读取实验设计、Feature表,并进行数据筛选和交叉筛选

# 读取实验设计和物种分类文件
# 之前有过下面一样的步骤,但确保正确在重新运行一次。

train = subset(metadata_f, Site_1 %in% c("S","E"))
train
idx3 = sample(1:386, size = 200, replace = F)
# 注size=67,66比较小
# 选择的样本标记为TRUE,未选择的为FALSE
idx3 = 1:386 %in% idx3
#train2 = metadata_f[idx3,]
train2 = train[idx3,] # train2是根据实验设计在train里面又筛选的作为训练集的数据。测试集是剩余部分。
train2
train2$Site_1 = as.factor(train2$Site_1)
summary(train2)
##    Filename          Group_day          Group_day1         Group_dayE       
##  Length:200         Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   Group_dayS           Group            StateSite         StateSite_1       
##  Length:200         Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##     State               Site           Site_1       Day           Day_num     
##  Length:200         Length:200         E: 79   Min.   : 0.00   Min.   : 0.00  
##  Class :character   Class :character   S:121   1st Qu.: 0.00   1st Qu.: 0.00  
##  Mode  :character   Mode  :character           Median : 3.00   Median : 3.00  
##                                                Mean   : 3.32   Mean   : 3.32  
##                                                3rd Qu.: 7.00   3rd Qu.: 7.00  
##                                                Max.   :10.00   Max.   :10.00  
##      Day1               DayE            Individual       
##  Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
## 
# 根据之前筛选情况读取物种分类水平数据。读取对应水平特征表
table =read.table(paste0("C:/16S/result/randomforestclassification/","tax_S_E_new_6Genus.txt"),header = T, row.names = 1)
# train2重命名,方便后面区分。
metadata_train1 = train2
summary(metadata_train1)
##    Filename          Group_day          Group_day1         Group_dayE       
##  Length:200         Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   Group_dayS           Group            StateSite         StateSite_1       
##  Length:200         Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##     State               Site           Site_1       Day           Day_num     
##  Length:200         Length:200         E: 79   Min.   : 0.00   Min.   : 0.00  
##  Class :character   Class :character   S:121   1st Qu.: 0.00   1st Qu.: 0.00  
##  Mode  :character   Mode  :character           Median : 3.00   Median : 3.00  
##                                                Mean   : 3.32   Mean   : 3.32  
##                                                3rd Qu.: 7.00   3rd Qu.: 7.00  
##                                                Max.   :10.00   Max.   :10.00  
##      Day1               DayE            Individual       
##  Length:200         Length:200         Length:200        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
## 
# 筛选OTU
idx = rownames(metadata_train1) %in% colnames(table)
metadata_train1 = metadata_train1[idx,]
otu_sub1 = table[, rownames(metadata_train1)] 
otu_sub1
dim(otu_sub1)
## [1]  92 200

选择最佳随机数(可选)

library(randomForest)
set.seed(77)
for (i in 0:9){
  set.seed(i)
  rf2 = randomForest(t(otu_sub1), metadata_train1$Site_1, importance=TRUE, proximity=TRUE, ntree = 5000)
  print(i)
  print(rf2)
}
## [1] 0
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 17.5%
## Confusion matrix:
##    E   S class.error
## E 50  29  0.36708861
## S  6 115  0.04958678
## [1] 1
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 17%
## Confusion matrix:
##    E   S class.error
## E 52  27  0.34177215
## S  7 114  0.05785124
## [1] 2
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 18%
## Confusion matrix:
##    E   S class.error
## E 52  27  0.34177215
## S  9 112  0.07438017
## [1] 3
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 16.5%
## Confusion matrix:
##    E   S class.error
## E 52  27  0.34177215
## S  6 115  0.04958678
## [1] 4
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 17%
## Confusion matrix:
##    E   S class.error
## E 53  26   0.3291139
## S  8 113   0.0661157
## [1] 5
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 17.5%
## Confusion matrix:
##    E   S class.error
## E 51  28  0.35443038
## S  7 114  0.05785124
## [1] 6
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 18%
## Confusion matrix:
##    E   S class.error
## E 52  27  0.34177215
## S  9 112  0.07438017
## [1] 7
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 17.5%
## Confusion matrix:
##    E   S class.error
## E 53  26  0.32911392
## S  9 112  0.07438017
## [1] 8
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 16.5%
## Confusion matrix:
##    E   S class.error
## E 52  27  0.34177215
## S  6 115  0.04958678
## [1] 9
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 17%
## Confusion matrix:
##    E   S class.error
## E 53  26   0.3291139
## S  8 113   0.0661157

随机森林分类

在确定的分类层级和最佳随机数下建模

library(randomForest)
set.seed(8)
rf3 = randomForest(t(otu_sub1), metadata_train1$Site_1, importance=TRUE, proximity=TRUE, ntree = 5000)
print(rf3)
## 
## Call:
##  randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000,      importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 5000
## No. of variables tried at each split: 9
## 
##         OOB estimate of  error rate: 16.5%
## Confusion matrix:
##    E   S class.error
## E 52  27  0.34177215
## S  6 115  0.04958678

交叉验证选择重要特征

set.seed(1) # 随机数据保证结果可重复,必须
# rfcv是随机森林交叉验证函数:Random Forest Cross Validation
result = rfcv(t(otu_sub1), metadata_train1$Site_1, cv.fold=5)
# 查看错误率表,31时错误率最低,为最佳模型
result$error.cv
##    92    46    23    12     6     3     1 
## 0.190 0.200 0.240 0.225 0.280 0.345 0.425
# 绘制验证结果 
with(result, plot(n.var, error.cv, log="x", type="o", lwd=2))

# 多次绘制
## 建立数据框保存多次结果
error.cv0 = data.frame(num = result$n.var, error.1 =  result$error.cv)
## 指定随机数循环5次
for (i in 1:(1+9)){
  print(i)
  set.seed(i)
  result= rfcv(t(otu_sub1), metadata_train1$Site_1, cv.fold=5) #  scale = "log", step = 0.9
  error.cv0 = cbind(error.cv0, result$error.cv)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
error.cv0 

绘制交叉验证曲线

# 提取x轴标签
n.var = error.cv0$num
# 提取y轴数据+标签
error.cv = error.cv0[,2:6]
colnames(error.cv) = paste('err',1:5,sep='.')
# 添加均值
err.mean = apply(error.cv,1,mean)
# 合并新的数据库,x+error+mean
allerr = data.frame(num=n.var,err.mean=err.mean,error.cv)
# number of otus selected 人为在图中观察的结果,30几乎为最低,且数量可接受
optimal = 20

# 图1:机器学习结果交叉验证图,选择Top features
# 图中 + 5条灰色拆线+1条黑色均值拆线+一条最优垂线+X轴对数变换
write.table(allerr, file = rfcv_filename_txt, sep = "\t", quote = F, row.names = T, col.names = T)

p = ggplot() + # 开始绘图
  geom_line(aes(x = allerr$num, y = allerr$err.1), colour = 'grey') + # 5次验证灰线 
  geom_line(aes(x = allerr$num, y = allerr$err.2), colour = 'grey') + 
  geom_line(aes(x = allerr$num, y = allerr$err.3), colour = 'grey') + 
  geom_line(aes(x = allerr$num, y = allerr$err.4), colour = 'grey') + 
  geom_line(aes(x = allerr$num, y = allerr$err.5), colour = 'grey') + 
  geom_line(aes(x = allerr$num, y = allerr$err.mean), colour = 'black') + # 均值黑线
  geom_vline(xintercept = optimal, colour='black', lwd=0.36, linetype="dashed") + # 最优垂线
  coord_trans(x = "log2") + # X轴对数变换和刻度
  scale_x_continuous(breaks = c(1, 2, 5, 10, 20, 30, 50, 100, 200)) + # , max(allerr$num)
  labs(title=paste('Training set (n = ', dim(t(otu_sub1))[1],')', sep = ''), 
       x='Number of families ', y='Cross-validation error rate') + 
  annotate("text", x = optimal, y = max(allerr$err.mean), label=paste("optimal = ", optimal, sep="")) + theme_bw()
p  

ggsave(p, file = rfcv_filename_pdf, width = 180, height = 120, unit = 'mm')

特征重要性可视化

## 预览和保存特征贡献度
imp= as.data.frame(rf3$importance)
imp = imp[order(imp$MeanDecreaseAccuracy, decreasing = T),]
# bbb = head(imp,n=optimal)
write.table(imp,file = importance_filename_txt,quote = F,sep = '\t', row.names = T, col.names = T)

# 读取所有feature贡献度
imp2 = read.table(importance_filename_txt, header=T, row.names= 1, sep="\t") 
# 分析选择top20分组效果最好,参数显示数量
imp2 = head(imp2, n = optimal)
imp2 = imp2[order(imp2$MeanDecreaseAccuracy, decreasing = F),]
# 简化全名,去掉界
imp2$Family = gsub("Bacteria\\|","",rownames(imp2))
# 添加门用于着色(删除竖线后面全部)
imp2$Phylum = gsub("\\|.*","",imp2$Family)

# 设置顺序
imp2$Family = factor(imp2$Family, levels = imp2$Family)

# 图2. 绘制物种类型种重要性柱状图
p = ggplot(imp2, aes(x = Family, y = MeanDecreaseAccuracy, fill = Phylum)) +   
  geom_bar(stat = "identity") + 
  coord_flip() + theme_bw()
p

ggsave(top_feautre_filename, p, width=200*2.5, height=59*2, unit='mm')
# 名称不一定唯一,需要手动修改

测试集独立验证

如果第一地点数据量足够大,可以取出1/2到1/3进行同一地点的独立验证。方法相同。 筛选测序集样品

# 测试集选择

test = train[!idx3,]
test
#test = metadata_f[!idx3,]
metadata_test = test
summary(metadata_test)
##    Filename          Group_day          Group_day1         Group_dayE       
##  Length:186         Length:186         Length:186         Length:186        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   Group_dayS           Group            StateSite         StateSite_1       
##  Length:186         Length:186         Length:186         Length:186        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##     State               Site              Site_1               Day        
##  Length:186         Length:186         Length:186         Min.   : 0.000  
##  Class :character   Class :character   Class :character   1st Qu.: 0.000  
##  Mode  :character   Mode  :character   Mode  :character   Median : 3.000  
##                                                           Mean   : 3.312  
##                                                           3rd Qu.: 7.000  
##                                                           Max.   :10.000  
##     Day_num           Day1               DayE            Individual       
##  Min.   : 0.000   Length:186         Length:186         Length:186        
##  1st Qu.: 0.000   Class :character   Class :character   Class :character  
##  Median : 3.000   Mode  :character   Mode  :character   Mode  :character  
##  Mean   : 3.312                                                           
##  3rd Qu.: 7.000                                                           
##  Max.   :10.000
idx4 = rownames(metadata_test) %in% colnames(table)
metadata_test = metadata_test[idx4,]
otu_sub2 = table[,rownames(metadata_test)]


# 转置,并添加分组信息
otutab_t2 = as.data.frame(t(otu_sub2))
otutab_t2$Site_1 = metadata_test[rownames(otutab_t2),]$Site_1

基于训练集随机森林模型验证

set.seed(315)
otutab.pred = predict(rf3, t(otu_sub2) )  
pre_tab = table(observed=otutab_t2[,"Site_1"],
                predicted=otutab.pred) 
pre_tab
##         predicted
## observed   E   S
##        E  54  22
##        S  10 100

可视化验证结果

# 整理样本原始分组和预测分类
predict = data.frame(group = otutab_t2[,"Site_1"], predicted=otutab.pred)

# 保存预测结果表
write.table("SampleID\t", file=prediction_filename,append = F, quote = F, eol = "", row.names = F, col.names = F)
write.table(predict, file = prediction_filename,append = T, quote = F, row.names = T, col.names = T, sep = "\t")
## Warning in write.table(predict, file = prediction_filename, append = T, :
## 给文件加列名
# 转换为数值可视化
# 预测准确标为1,错误标为0
predict$result = ifelse(predict$group == predict$predicted, 1, 0)
# IND=1, TEJ=2
predict$predict = ifelse(predict$predicted == "E", 1, 2)
# Set sample number in each row
column = 10
low_to_high_colors <- colorRampPalette(c("#B1DDE1", "#F8C7BD"))(100)

# Draw prediction result
AA1 = predict[predict$group=="S",]$predict
length(AA1)
## [1] 110
row = round(length(AA1)/column + 0.5)
row
## [1] 12
i = column * row - length(AA1)
AA1 = c(AA1, rep(NA, i))
matrix1 = matrix(AA1, nrow = row, ncol = column, byrow = T)
pheatmap(matrix1, cluster_rows = F, cluster_cols = F, cellwidth = 15, cellheight = 12, color = low_to_high_colors, border_color = "#EFEFEF")

pheatmap(matrix1, cluster_rows = F, cluster_cols = F, cellwidth = 18, cellheight = 18,filename = testfilename1, color = low_to_high_colors, border_color = "#EFEFEF")

# Draw another prediction result
BB1 = predict[predict$group=="E",]$predict
length(BB1)
## [1] 76
row = round(length(BB1)/column + 0.5)
i = column * row - length(BB1)
BB1 = c(BB1, rep(NA, i))
matrix2 = matrix(BB1, nrow = row, ncol = column, byrow = T)
pheatmap(matrix2, cluster_rows = F, cluster_cols = F, cellwidth = 15, cellheight = 12, color = low_to_high_colors, border_color = "#EFEFEF")
# 保存图片
pheatmap(matrix2, cluster_rows = F, cluster_cols = F, cellwidth = 18, cellheight = 18, filename = testfilename2, color = low_to_high_colors, border_color = "#EFEFEF")


## 将生成的文件放入指定文件夹
source_dir <- "C:/16S/result/randomforestclassification"
target_dir <- result_file_names

if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}

# 选取含特定字符的文件(例如后缀为 csv 的所有文件)
tobeCopy1 <- list.files(source_dir, pattern = "*_rf.txt", full.names = TRUE)
tobeCopy2 <- list.files(source_dir, pattern = "*_rf.pdf", full.names = TRUE)

# 复制选中文件到目标文件夹
sapply(tobeCopy1, function(x) { file.copy(x, target_dir, overwrite = TRUE) })
##        C:/16S/result/randomforestclassification/new_importance_S_E_rf.txt 
##                                                                      TRUE 
## C:/16S/result/randomforestclassification/new_prediction_binary_S_E_rf.txt 
##                                                                      TRUE 
##              C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.txt 
##                                                                      TRUE 
##          C:/16S/result/randomforestclassification/sub_metadata_S_E_rf.txt 
##                                                                      TRUE
sapply(tobeCopy2, function(x) { file.copy(x, target_dir, overwrite = TRUE) })
##         C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.pdf 
##                                                                 TRUE 
##           C:/16S/result/randomforestclassification/new_test_E_rf.pdf 
##                                                                 TRUE 
##           C:/16S/result/randomforestclassification/new_test_S_rf.pdf 
##                                                                 TRUE 
## C:/16S/result/randomforestclassification/new_top_feautres_S_E_rf.pdf 
##                                                                 TRUE
# 删除源文件夹中的已复制文件
file.remove(tobeCopy1)
## [1] TRUE TRUE TRUE TRUE
file.remove(tobeCopy2)
## [1] TRUE TRUE TRUE TRUE
# 确认文件已经复制和删除
print("Files copied and deleted successfully.")
## [1] "Files copied and deleted successfully."

使用此脚本,请引用下文:

If used this script, please cited:

Yong-Xin Liu, Lei Chen, Tengfei Ma, Xiaofang Li, Maosheng Zheng, Xin Zhou, Liang Chen, Xubo Qian, Jiao Xi, Hongye Lu, Huiluo Cao, Xiaoya Ma, Bian Bian, Pengfan Zhang, Jiqiu Wu, Ren-You Gan, Baolei Jia, Linyang Sun, Zhicheng Ju, Yunyun Gao, Tao Wen, Tong Chen. 2023. EasyAmplicon: An easy-to-use, open-source, reproducible, and community-based pipeline for amplicon data analysis in microbiome research. iMeta 2: e83. https://doi.org/10.1002/imt2.83

Copyright 2016-2023 Yong-Xin Liu , Tao Wen , Tong Chen